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EXECUTIVE  SUMMARY 


A  PC  based  three  dimensional  computational  procedure  has  been  developed.  With  this 
procedure,  the  engineer  will  be  able  to  check  the  design  of  funnel  and  gate  systems  using  the  site 
specific  properties  and  conditions. 

The  name  for  the  PC  Smooth  Particle  Hydrodynamics  (SPH)  funnel  and 
gate/groundwater  flow  computer  program  is  GROWFLOW.  The  code  is  presently  fully 
operational  on  a  486  DX2  PC  running  under  the  LINUX  operating  system  and  has  the  following 
capabilities: 

1.  Fully  3D 

2.  Handles  both  saturated  and  partially  saturated  flow  conditions 

3.  Saturated  hydraulic  conductivity  can  vary  throughout  the  flow  domain  (can  map  very 
complex  in  situ  conditions) 

4.  Saturated  hydraulic  conductivity  can  be  orthotopic 

5.  In-flow  conditions  can  vary  with  time 

6.  Handle  many  sources  simultaneously 

7.  Arbitrarily  complex  saturation/head  and  saturation/hydraulic  conductivity 
relationships  using  look-up  tables 

8.  Features  exterior  and  interior  flow  control  panels  to  direct  and  contain  flow  (build  any 
complex  shape  desired) 

9.  No  "grid"  required  (only  have  to  place  the  integration  points  so  complex  shapes  are 
easy  to  model) 
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SECTION  I 


INTRODUCTION 


OBJECTIVE 

The  overall  objective  of  this  program  is  to  develop  a  three-dimensional  analysis  tool  for 
the  engineer  engaged  in  the  design  and  evaluation  of  funnel  and  gate  systems. 

BACKGROUND 

Funnel-and-gate  systems  are  one  approach  to  the  remediation  of  the  groundwater  which 
has  been  contaminated.  These  funnels  (an  impervious  wall)  are  used  to  direct  the  contaminated 
groundwater  flow  into  a  gate  region  where  remediation  can  take  place.  The  system  serves  two 
basic  purposes.  First  the  system  directs  the  flow  and  thereby  contains  the  contaminant  plume 
and  secondly  the  contaminant  in  the  groundwater  is  at  higher  concentrations. 

SCOPE 

The  scope  of  this  study  is  to  develop  a  PC-based  three-dimensional  computational 
procedure  for  the  analysis  and  design  of  funnel  and  gate  systems.  The  approach  used  is  a 
modified  form  of  Smooth  Particle  Hydrodynamics  (SPH).  The  resulting  computer  program  is 
named  GROWFLOW. 

The  GROWFLOW  code  is  extremely  flexible.  Complex  in  situ  conditions  and 
geometries  are  easily  accommodated.  Another  important  feature  of  the  code  is  that  it  can  model 
highly  nonlinear  systems  very  efficiently.  Convergence  is  guaranteed.  This  code  operates  on  a 
totally  different  time  regune  than  a  typical  hydrocode.  The  stable  time  step  for  a  hydrocode 
might  be  on  the  order  of  tens  to  himdreds  of  microseconds.  The  stable  time  step  for 
GROWFLOW  depends  on  the  hydraulic  conductivity  of  the  media.  The  following  report 
summarizes  the  theoretical  base  for  the  procedure. 
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SECTION  II 


GROWFLOW  SOLUTION  ALGORITHMS 

In  the  following  sections,  the  solution  algorithms  for  the  GROWFLOW  code  will  be 
described.  First,  the  explicit  time-integration  solution  method  will  be  discussed,  and  a  flowchart 
of  the  time  integration  loop  will  be  presented.  Next,  the  SPH  (Smooth  Particle  Hydrodynamics) 
technique  that  the  analysis  method  is  based  on  will  be  presented.  Then,  the  partial  differential 
equations  (conservation  of  mass  and  momentum)  will  be  derived  and  will  be  cast  into  the  SPH 
formalism.  Since  the  explicit  time  integration  method  is  conditionally  stable,  a  formulation  for 
the  propagation  speed  in  the  flow  field  will  be  presented,  along  with  a  derivation  of  the  stability 
limit  for  the  method.  A  method  for  controlling  external  and  internal  flow  patterns  (flow  control 
panels)  will  be  explained.  The  methods  for  treating  source  regions  will  be  described.  Finally, 
the  method  used  for  describing  the  constitutive  relationships  for  the  flow  media  will  be 
presented. 

TIME  INTEGRATION  LOOP 

The  GROWFLOW  analysis  code  uses  an  explicit  time-integration  solution  method  that 
will  be  described  in  this  section.  The  basis  for  the  time-integration  method  will  be  given,  and  a 
flowchart  of  the  time  integration  loop  will  be  presented. 

An  explicit  time-integration  solution  method  was  chosen  for  the  GROWFLOW  code. 
There  are  a  few  reasons  why  this  is  a  superior  choice  for  the  GROWFLOW  code.  First,  the 
explicit  time-integration  method  is  the  method  of  choice  for  SPH  based  techniques.  This  is 
necessary  because  of  the  nature  of  the  SPH  method.  Coimectivity  between  integration  points  is 
loosely  defined,  making  the  assembly  of  the  matrices  necessary  for  implicit  time-integration 
techniques  difficult.  Next,  the  behavior  of  the  groimdwater  flow  problem  is,  in  general,  highly 
nonlinear.  Because  the  explicit  time  integration  method  uses  small  time  steps,  simulation  of  non¬ 
linear  effects  does  not  require  any  additional  treatment  beyond  the  normal  time-stepping  (i.e. 
subiterations  are  not  required).  Finally,  the  explicit  time  integration  method  eliminates  the  need 
for  the  solution  of  large  sets  of  simultaneous  equations.  The  system  of  partial  differential 
equations  is  solved  sequentially  instead  of  simultaneously.  This  can  be  a  big  advantage  for 
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solving  large  three-dimensional  problems.  The  memory  requirement  for  the  solution  is  greatly 
reduced.  This  allows  for  solution  on  a  PC  class  machine,  versus  a  large  workstation  or 
supercomputer.  Details  on  the  explicit  time  integration  method  follow. 

The  explicit  time-integration  method  is  a  time-stepping  method.  Time  steps  are  chosen 
so  that  variables  change  by  very  small  amounts  each  time  step.  Variables  (functions)  are 
updated  at 

each  time  step  using  the  following  formula: 


F 


F''  + 


(1) 


where, 

F'  =  Function  value  at  present  time  step 
F'**  =  Function  value  at  previous  time  step 
A  t  =  Time  step  increment 

The  selection  of  a  stable  time  step  is  critical  to  the  method  and  will  be  explained  in  a  later 
section.  For  each  small  time  increment,  the  solution  is  performed  as  shown  in  Figure  1 . 

SPH  SOLUTION  ALGORITHM 

The  GROWFLOW  analysis  code  uses  a  unique  Lagrangian  solution  technique.  The 
techmque  is  based  on  SPH  techniques.  For  excellent  reviews  of  the  SPH  technique  and  its 
background  see  Monaghan,  1992  and  Benz,  1989.  A  discussion  of  the  advantages  of  the  SPH 
method  for  groundwater  flow  simulation  will  be  presented  first.  A  brief  description  of  the 
original  SPH  formulation  ^11  then  be  given.  The  SPH  technique  that  is  used  in  GROWFLOW  is 
completely  new  and  original.  A  description  of  this  new  SPH  method,  called  “Gradient  SPH,” 
will  be  given. 

The  SPH  method  has  some  very  important  advantages  for  groundwater  flow  simulation. 
Probably  the  single  most  important  advantage  is  that  it  is  a  fully  Lagrangian  method.  Almost  all 
other  groundwater  flow  analysis  codes  are  Eulerian.  The  important  advantage  for  Lagrangian 
codes  is  that  they  track  material  points.  This  eliminates  the  expensive  “advection”  step  that 
Eulerian  codes  must  perform.  It  also  eliminates  the  “dispersion”  effect  often  seen  in  Eulerian 
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codes.  Sharp,  clean  material  interfaces  are  preserved. 

The  SPH  method  involves  the  use  of  arbitrarily  arranged  integration  points  that  interact 
with  each  other  through  an  interpolation  function  or  interpolation  kernel.  The  integration  points 
are  dispersed  throughout  the  body  to  provide  a  continuum  description  of  the  body.  The 
technique  is  strictly  Lagrangian,  since  integration  points  represent  material  points.  To  aid  in  the 
discussion,  see  Figure  2.  Integration  point  2  is  within  the  influence  zone  of  integration  point  1 . 
Therefore,  they  can  interact.  Their  interaction  is  controlled  by  their  separation  distance,  R,  and 
the  smoothing  length  (typically  denoted  by  h).  The  value  of  the  interpolation  function,  W,  is  a 
function  of  R  and  the  smoothing  length.  There  are  only  two  integration  points  shown  in  the 
figure,  but  an  integration  point  will  typically  have  many  different  neighbors  within  its  influence 
zone.  Each  neighbor  within  the  influence  zone  contributes  to  the  response  of  the  integration 
point  being  considered.  By  summing  up  all  of  the  interactions,  the  behavior  of  the  whole  body 
can  be  approximated.  This  can  be  expressed  mathematically  with  a  summation  equation: 

Fi  =  Fj  —  W  (2) 

where, 

F  =  Any  scalar  or  vector  valued  function 
i  j  =  Summation  indices 

N  =  Total  number  of  integration  points  that  are  within  influence  zone  of  integration 
point  i 
m  =  Mass 
p  =  Density 

W  =  Kernel  (interpolation  polynomial) 

The  kernel  can  have  many  different  forms.  Refer  to  Monaghan,  1992  or  Benz,  1989  for  a 
discussion  on  the  rules  for  choosing  and  constructing  kernels.  The  kernel  that  was  chosen  for  the 
GROWFLOW  formulation  is  an  exponential  kernel  (Benz,  1989).  It  has  the  following  form  for 
the  three-dimensional  case: 


where, 

h  Smoothing  length 

V  R/h 

A  completely  original  form  of  the  SPH  technique  was  developed  in  this  effort.  Early 
versions  of  GROWFLOW  used  traditional  SPH  techniques.  Results  were  not  satisfactory  in 
large  flow  regions,  so  a  new  method  was  developed  that  would  be  extremely  accurate  for  any 
general  groundwater  flow  regime.  A  new  method  called  “Gradient  SPH”  was  developed  for  the 
GROWFLOW  code.  As  the  name  implies,  the  method  focuses  in  on  the  accurate  solution  of 
equations  involving  gradients  and  differentials.  The  partial  differential  equations  of  motion  used 
to  model  groundwater  flow  only  involve  equations  of  this  type.  Consider  the  set  of  SPH 
equations  shown  above  to  understand  the  need  for  such  a  formulation. 

The  need  for  a  gradient-based  SPH  method  can  be  explained  by  considering  the  two  SPH 
Equations  (2)  and  (3).  The  SPH  method  is  a  volume  based  method.  The  mass-over-density  term 
in  Equation  (2)  is  called  the  number  density  function.  It  is  the  way  weighting  is  done  in  SPH. 
The  number  density  function  produces  a  volume  term.  Consideration  of  the  kernel  function  in 
Equation  (3)  shows  that  it  produces  a  one  over  volume  term.  This  volume  weighting  procedure 
works  extremely  well  for  areas  inside  of  the  body  (away  from  boundaries).  Areas  near 
boundaries  present  a  problem.  The  description  of  the  volume  is  incomplete,  since  there  are  no 
integration  points  to  sum  over  outside  of  the  boundary.  Groundwater  flow  problems  can  be 
dominated  by  boundary  effects.  Clearly  a  method  was  needed  which  is  highly  accmate  in  the 
interior  and  near  the  exterior  of  the  solution  domain.  **Gradient  SPH”  was  developed  in  response 
to  that  need. 

The  idea  behind  Gradient  SPH  is  simple,  yet  powerful.  Consider  the  first  derivative  of 
the  SPH  function  in  Equation  (2)  with  respect  to  position: 

F;  =  ZjL;  Fj  ^  VW  (4) 

Replace  the  number  density  function  in  Equation  (4)  with  the  original  volume  and  add  a 
normalization  function: 
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f;  •  s”,  Fj  V.,  vw 


(5) 


where, 

C,  =  Normalization  function 
Vo  =  Original  volume 

The  normalization  function  is  defined  to  preserve  the  gradient,  so  that  the  gradient  will  be 
accurately  reproduced  no  matter  how  many  neighbors  the  integration  point  has  (eliminate 
boundary  effects).  A  normalization  function  that  produces  this  effect  is; 

=  [sj'.,  V.,  (Xi-Xj)>vw]  '  (6) 

where, 

X  =  Cartesian  coordinates  (x  =  x,  y,  z) 

The  set  of  Equations  (5)  and  (6)  are  the  set  of  SPH  equations  that  make  up  the  Gradient  SPH 
method.  They  are  the  basic  SPH  equations  used  in  GROWFLOW.  An  exponential  kernel  is 
used  in  GROWFLOW.  The  gradient  of  the  kernel  function  is  defined  as: 

Vlf  =  ^  ^  exjl'j}  (7) 

h  ox 

In  order  to  keep  the  method  symmetric,  the  smoothing  length  is  defined  as: 

h  =  ^  (h  +  hj)  (8) 

A  nonsymmetric  interpolation  kernel  can  be  used  in  GROWFLOW.  This  is  extremely 
useful  in  flow  fields  where  the  flow  is  compressing  in  one  direction  and  expanding  in  another.  In 
areas  where  the  flow  is  compressing,  greater  accuracy  can  be  achieved  by  allowing  the 
smoothing  length  to  shrink.  In  areas  where  the  flow  is  expanding,  greater  accuracy  can  be 
achieved  by  allowing  the  smoothing  length  to  increase.  This  is  equivalent  to  keeping  the  number 
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zone  of  influence  that  is  an  ellipsoid.  The  time  rate  of  change  of  the  smoothing  length  in  a 
particular  direction  follows  the  formulation  given  by  Benz,  1989: 


dt 


(9) 


where, 

=  Velocity  in  a  given  direction 
GOVERNING  PARTIAL  DIFFERENTIAL  EQUATIONS 

The  governing  partial  differential  equations  for  groundwater  flow  in  a  Lagrangian  frame 
of  reference  will  be  derived  and  will  be  cast  into  the  SPH  formalism  in  this  section.  The  two 
partial  differential  equations  of  interest  are  the  conservation  of  mass  equation  (continuity 
equation)  and  the  conservation  of  momentum  equation  (Darcy’s  Law).  The  continuity  equation 
will  be  derived  first  and  cast  into  the  SPH  formalism.  The  continuity  equation  will  be  derived 
for  two  different  flow  conditions:  partially  saturated  and  fully  saturated.  Then  the  conservation 
of  momentum  equation  will  be  presented  and  cast  into  the  SPH  formalism. 

The  conservation  of  mass  (continuity  equation)  will  be  derived  and  cast  into  the  SPH 
formalism.  Begin  by  defining  some  terms  that  will  be  used  in  the  derivation: 
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P 

S 


=  moisture  content  = 

=  fluid  density  =  — 
V  w 


K 

Vr 


=  Degree  of  Saturation  = 


n 


where, 

V„  =  Volume  of  fluid 

V-r  =  Total  volume  (volume  of  soil  plus  volume  of  fluid) 
m^^  =  Mass  of  fluid 
n  =  Porosity 
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To  proceed  with  the  derivation  of  the  conservation  of  mass,  consider  the  differential  volume 
shown  in  Figure  3.  The  mass  of  water  in  this  differential  element  is  defined  as: 


niy,  =  p  Odxdydz 


(10) 


Take  the  derivative  of  Equation  (10)  with  respect  to  time  to  get  the  time  rate  of  change  of  the 
mass  of  fluid  inside  of  the  differential  element: 


dmv> 

dt 


d(pe) 

dt 


dxdydz 


(11) 


The  mass  flow  rate  of  fluid  through  the  boundary  of  the  differential  element  in  a  given  direction 
is  defined  as: 


S(P&J 

dx 


dxA 


d(p&J 

dx 


dxdydz 


(12) 


Or,  the  sum  of  the  mass  flow  rates  through  the  boundaries  of  the  differential  element  is: 

W •(p&^)  dx  dy  dz  (13) 

The  conservation  of  mass  states  that  the  negative  of  the  time  rate  of  change  of  the  fluid  mass  in 
the  differential  element  must  equal  the  sum  of  the  fluid  mass  flow  rates  through  the  boundaries  of 
the  differential  element.  Equation  (11)  must  equal  Equation  (13): 


dt 

Simplify  Equation  (14): 


d(pe) 

dt 


V*(p3J 


(15) 
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(15) 


d{pe) 

dt 


Assume  that  the  fluid  density  does  not  vary  with  time  or  position  (the  fluid  in  incompressible). 
The  fluid  density  term  can  be  moved  outside  of  the  derivatives  in  the  right-hand  side  and  left- 
hand  side  of  Equation  (15): 


(16) 


The  fluid  density  can  now  be  eliminated  from  both  sides  of  Equation  (16): 


-—  = 
dt 

Recall  the  following  identity: 


e  =  nS 


(17) 


(18) 


Substitution  of  Equation  (18)  into  Equation  (17)  yields  the  final  form  of  the  conservation  of  mass 
(continuity)  equation: 


d(nS) 

dt 


V.9, 


(19) 


The  conservation  of  mass  equation  can  be  used  to  calculate  the  rate  at  which  either  the 
degree  of  saturation,  S,  or  the  porosity,  n,  changes  with  time.  There  are  two  flow  cases  to  be 
considered  here.  For  the  partially  saturated  flow  condition,  the  porosity  will  not  change  with 
time.  For  the  fully  saturated  flow  condition,  the  degree  of  saturation  will  not  change  with  time. 
If  these  conditions  are  imposed  on  Equation  (19),  the  following  set  of  two  equations  are  derived: 


dS 

1 

dt 

—  v.d 

n 

Partially  Saturated 

(20) 

dn 

1 

Fully  Saturated 

dt 

— 

S'  * 

(21) 

9 


r  (1  ^ 

f 

< 

• 

1 

50 

-  ^x* 

v| 

- 

\n  V 

KnJJ 

Convert  Equation  (22)  to  the  Gradient  SPH  form  in  Equation  (5): 


dSi 

dt 


Vtj 


•vw 


njJ 


(22) 


(23) 


Simplify  Equation  (23)  to  get  the  equation  used  to  calculate  the  rate  of  change  of  the  degree  of 
saturation  for  partially  saturated  flow  regions: 


5*?  Vt  l  \ 

Partially  Saturated 


(24) 


Next,  consider  the  fully  saturated  condition.  Manipulate  Equation  (21)  using  the  chain 

rule: 


dn 

dt 


-  &A  -  &• 

V  - 

vs  ' 

V  \SJ) 

(25) 


Convert  Equation  (25)  to  the  Gradient  SPH  form  in  Equation  (5): 


drii 


-c  SjL/  Vtj 


&x  ^ 

Xi 

Sj'SjJ 


•vw 


(26) 


Simplify  Equation  (26)  to  get  the  equation  used  to  calculate  the  rate  of  change  of  porosity  for 
fully  saturated  flow  regions: 


10 


=  4',  SjL/  (>9x,. Fully  Saturated 


(27) 


Equations  (24)  and  (27)  are  the  conservation  of  mass  (continuity)  equations  used  in  the 
GROWFLOW  method.  The  conservation  of  momentum  equation  will  be  discussed  next. 

The  conservation  of  momentum  equation  will  be  presented  and  cast  into  the  SPH 
formalism.  The  conservation  of  momentum  equation  used  is  a  modified  form  of  Darcy’s  Law, 
suggested  by  Huyakom  et  al.,  1986: 

S,  =  -kK.S/(<p  +  Z)  (28) 

where, 

k  =  Relative  permeability  with  respect  to  the  water  phase  (function  of  S) 

K  =  Saturated  hydraulic  conductivity  tensor 
^  =  Pressure  head 
Z  =  Elevation  head 

The  conservation  of  momentum  equation  can  be  cast  into  the  SPH  formalism  as  follows. 
Rearrange  Equation  (28)  using  the  chain  rule: 

<9^=  -k[V»(y/K)  -  [{N*K]  (29) 

where, 

xj/  -  Total  head  =  ^  +  Z 

Convert  Equation  (29)  to  the  Gradient  SPH  form  in  Equation  (5): 

=  -k^  C  Vtj  (KjV^j  -  (30) 

Simplify  Equation  (30)  to  get  the  form  of  the  conservation  of  momentum  equation  used  in  the 
GROWFLOW  method: 
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=  ki  Vt,  (w,  -  ¥^j)  Kj*VW 


(31) 


The  set  of  three  Equations  (24),  (27),  and  (31)  form  the  set  of  equations  that  are  used  to 
solve  the  governing  partial  differential  equations  in  the  GROWFLOW  method. 

PROPAGATION  SPEED  AND  STABLE  TIME  STEP 

A  formulation  for  the  propagation  speed  of  a  disturbance  in  a  groundwater  flow  field  will 
be  presented  in  this  section.  The  formulation  will  be  cast  into  the  SPH  formalism  for  use  in 
calculating  the  stability  limit  of  the  solution.  The  stability  limit  will  be  discussed,  and  the 
selection  criteria  for  time  step  size  will  be  given. 

The  following  is  a  presentation  of  the  formulation  for  the  propagation  speed  of  a 
disturbance  in  a  groimdwater  flow  field.  The  partial  differential  equation  is  derived  in  a 
Lagrangian  frame  of  reference.  Start  with  the  governing  partial  differential  equation. 

The  governing  partial  differential  equation  is  a  modified  form  of  Richard’s  equation 
(Therrien,  1992).  This  equation  is  obtained  by  substituting  Equation  (28),  the  conservation  of 
momentum  equation,  into  Equation  (17),  the  conservation  of  mass  equation: 


8t 


(32) 


Write  Equation  (32)  in  indicial  notation  and  split  the  differential  of  the  total  head: 


de 


dt  dXa 


k  K, 


dy/  66 


ap 


dO  dxp) 


a,p  =  x,y,z 


(33) 


Make  for  following  definition: 


dw 


(34) 


Substitute  Equation  (34)  into  Equation  (33): 
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SXa\  dXfi, 


Split  the  partial  differential  of  moisture  content  with  respect  to  time: 


dxa  ^  [  SO 

dt  dXa  dXa\“^  dxp. 


The  speed  of  propagation  in  a  particular  direction  is  defined  as: 


Substitute  the  definition  in  Equation  (37)  into  Equation  (38): 


r  —  -  —\ 

dXa  dxa\“^dxp 


Expand  Equation  (38)  using  the  chain  rule: 


^  dXp  dXa  dXp  dXa  dXp 


Eliminate  the  first  partial  differential  of  moisture  content  from  Equation  (39): 


C  =  ^  +  n  ^ 
^  dXa  “^dXa 


Insert  Equation  (34)  into  Equation  (40): 


Insert  Equation  (18)  into  Equation  (41),  and  note  that  head  is  not  a  function  of  porosity: 
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(42) 


dxa^n  dS)  n  dS  dxa 


Rewrite  Equation  (42)  in  tensor  notation: 


C  =  to; 


(43) 


Equation  (43)  is  the  partial  differential  equation  for  speed  of  propagation.  The  derivative  of  head 
with  respect  to  degree  of  saturation  is  a  constitutive  relationship  "nd  will  be  discussed  in  a  later 
section.  The  last  step  is  to  cast  it  into  the  SPH  formalism. 

The  speed  of  propagation  equation  can  be  cast  into  the  SPH  formalism  as  follows. 
Rearrange  Equation  (43)  using  the  chain  rule: 


C  =  AT.VI-  ^1  + 
n  dS 


k  dy/  „  „ 

— ^  v#a:  +  V* 

n  dS 


kKS 


dyf 

Is 


,  k  dw 
-  nS^*\-K^ 
clS 


Convert  Equation  (44)  to  the  Gradient  SPH  form  in  Equation  (5): 


(44) 


Q  =  sjiy  Vtj 


Kr 


dS 


.  ^  ^^~nSdS)^ 


dy/ 


\dS 


•  VJV 
(45) 


Equation  (45)  is  the  final  form  of  the  equation  used  to  calculate  speed  of  propagation.  The  speed 
of  propagation  can  be  used  to  calculate  the  stability  limit  of  the  calculation. 

The  explicit  time  integration  used  in  GROWFLOW  is  “conditionally  stable”.  A  time  step 
increment  must  be  used  that  is  smaller  than  the  stability  limit  for  the  solution.  If  the  time  step 
increment  is  too  big,  the  solution  will  become  unstable  and  numerical  oscillations  will  grow  very 
quickly.  Fortunately,  there  are  rules  for  estimating  the  stability  limit  for  a  calculation. 

The  stability  limit  can  be  estimated  for  a  particular  solution.  The  stability  limit  for  an 
explicit  time  integration  solution  is  generally  the  time  it  takes  for  a  disturbance  to  travel  a 
“characteristic  length”.  A  characteristic  length  for  a  finite  element  solution  is  the  length  of  the 
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The  stability  limit  can  be  estimated  for  a  particular  solution.  The  stability  limit  for  an 
explicit  time  integration  solution  is  generally  the  time  it  takes  for  a  disturbance  to  travel  a 
“characteristic  length”.  A  characteristic  length  for  a  finite-element  solution  is  the  length  of  the 
side  of  an  element.  In  SPH  solutions,  it  is  the  smoothing  length,  h.  No  time  step  increment  can 
be  larger  than  the  time  it  takes  a  disturbance  to  travel  a  smoothing  length.  This  can  be  stated 
mathematically: 

h 

^  ^  Q  (46) 

The  procedure  used  in  the  GROWFLOW  code  is  to  find  the  smallest  combination  of  smoothing 
length  over  speed  of  propagation,  and  then  multiply  this  by  a  safety  factor.  The  time  step 
increment  is  calculated  using  the  following  equation: 


where, 

X  =  Factor  of  safety  (0.5  is  used  in  GROWFLOW) 

A  new  time  step  increment  is  calculated  every  time  step  using  Equation  (47). 
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SECTION  III 


FLOW  CONTROL  PANELS 

Flow  control  panels  will  be  discussed  in  this  section.  Flow  control  panels  are  used  in 
GROWFLOW  to  limit  flow  at  the  external  boundaries  and  to  create  internal  boundaries  in  the 
solution.  The  characteristics  of  the  flow  control  panels,  how  they  operate,  and  were  they  may  be 
used  will  be  explained  in  the  following  discussion. 

Flow  control  panels  have  the  following  characteristics.  They  are  quadrilateral  shapes  that 
are  described  by  four  comers.  They  are  assumed  to  have  no  thickness.  Any  arbitrary  shape  can 
be  assembled  with  the  panels  (similar  to  using  finite  element  shell  or  plate  elements  to  describe  a 
shape).  The  panels  are  used  to  control  flow  through  the  area  within  the  limits  of  the  panel  edges. 

The  flow  control  panels  operate  by  preventing  flow  normal  to  the  surface  of  the  panel 
from  occurring.  They  form  an  impermeable  boundary.  They  can  be  used  on  the  external 
boundaries  to  prevent  fluid  from  flowing  out.  This  is  equivalent  to  putting  sides  on  a  tank.  They 
can  be  used  in  the  interior  to  prevent  flow  through  an  area.  Flow  parallel  to  the  panel  is  not 
restricted. 

The  algorithm  for  the  panels  works  in  the  following  way.  Fluid  is  not  allowed  to 
penetrate  the  panel.  All  of  the  fluid  particles  that  are  within  a  small  capture  zone  near  the  panel 
are  checked  for  penetration  at  each  time  step.  If  the  motion  of  the  particle  would  result  in 
penetration  of  the  panel,  the  velocity  component  normal  to  the  panel  is  set  to  zero.  This  prevents 
penetration  while  still  allowing  motion  parallel  to  the  panel.  Also,  the  penetration  point  must  be 
within  the  boundaries  of  the  edges  of  the  panel. 

The  flow  control  panels  can  be  used  on  the  exterior  or  in  the  interior  of  the  problem.  A 
separate  type  of  panel  is  used  for  each  situation.  The  interior  panel  differs  from  the  exterior 
panel  in  one  very  important  way.  Interior  panels  do  not  allow  interactions  to  occur  between 
particles  on  opposite  sides  of  the  panel.  Exterior  panels  do  not  check  for  this  condition  and 
should  not  be  used  in  the  interior  of  the  problem.  The  need  to  check  for  interactions  and  limit 
interactions  is  the  result  of  the  nature  of  SPH.  Recall  that  SPH  integration  points  (particles) 
interact  with  all  of  the  neighbors  within  an  influence  zone  (two  times  the  smoothing  length). 
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This  occurs  through  a  summation  (integration)  process.  An  interaction  between  two  neighbors 
separated  by  a  flow  control  panel  must  be  prevented  to  model  physical  reality. 

Flow  control  panels  are  an  excellent  tool  for  modeling  the  physical  impermeable 
boimdaries  of  a  problem.  They  can  also  be  used  to  precisely  model  internal  impermeable 
boimdaries  such  as  walls. 
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SECTION  IV 


SOURCE  MODELS 

The  two  source  models  available  in  GROWFLOW  will  be  described  in  this  section.  They 
consist  of  a  constant  head  source  model  and  a  constant  velocity  source  model.  The  common 
theory  behind,  and  operation  of,  the  source  models  will  be  discussed  first,  then  each  model  will 
be  described. 

A  special  procedure  must  be  used  for  the  source  models  in  GROWFLOW.  This 
procedure  consists  of  defining  the  integration  points  (particles)  that  make  up  the  source,  defining 
the  source  panel,  and  injecting  particles  from  the  source  into  the  flow  field.  Recall  that  the  SPH 
method  is  a  Lagrangian  method.  Material  points  are  tracked  in  the  solution.  All  of  the  material 
in  the  solution  must  be  defined  at  the  beginning  of  the  solution.  That  means  that  the  material  that 
flows  in  from  the  source  regions  must  be  accounted  for  (defined)  at  the  beginning  of  the  solution 
(see  Figure  4).  A  set  of  particles  are  defined  as  source  particles  (shown  in  red).  The  source 
particles  move  through  the  source  region  perpendicular  to  the  source  panel.  The  source  panel  is  a 
quadrilateral  defined  by  four  points,  just  like  flow  control  panels.  As  the  particle  approaches  the 
source  panel,  a  check  is  made  to  see  if  the  source  particle  penetrates  the  source  panel.  If 
penetration  occurs,  injection  of  the  source  particle  into  the  flow  field  occurs.  The  particle  is 
converted  to  a  normal  flow  field  particle.  Particles  in  the  source  region  can  only  interact  with 
regular  flow  field  particles  through  the  source  panel.  This  was  done  to  limit  interactions  to  just 
the  area  around  the  source  panel.  This  is  of  practical  importance  for  modeling  interior  sources  in 
the  flow  field. 

Note  that  the  source  models  in  GROWFLOW  only  allow  one-way  flow.  That  is,  only 
source  peirticles  can  pass  through  the  source  panel.  Regular  flow  field  particles  cannot  move  into 
the  source  region.  This  simply  means  that  the  source  models  should  not  be  used  to  model  sinks. 
A  better  way  to  model  sinks  in  GROWFLOW  is  to  simply  provide  a  volume  for  fluid  to  flow 
into.  There  are  two  types  of  source  models  in  GROWFLOW. 

A  constant  head  source  model  is  available  in  GROWFLOW.  Inflow  velocity  is  adjusted 
so  that  a  constant,  predefined  head  is  maintained  in  the  source  region. 
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A  constant  velocity  source  model  is  available  in  GROWFLOW.  The  head  in  the  source 
region  is  adjusted  so  that  a  constant,  predefined  inflow  rate  is  maintained. 
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SECTION  V 


CONSTITUTIVE  MODELS 

The  constitutive  models  used  in  GROWFLOW  will  be  described  in  this  section.  There 
are  two  constitutive  relationships  that  must  be  defined  for  a  groundwater  flow  analysis:  the 
relationship  between  degree  of  saturation  and  pressure  head,  and  the  relationship  between  degree 
of  saturation  and  relative  permeability.  These  two  relationships  will  be  discussed,  then  the 
method  for  describing  these  functions  in  GROWFLOW  will  be  explained. 

There  are  many  different  functional  relationships  that  have  been  used  to  describe  the 
relationship  between  degree  of  saturation,  S,  and  pressure  head,  n,  and  the  relationship  between 
degree  of  saturation,  S,  and  relative  permeability,  k.  As  Therrien  (1992)  points  out,  any 
physically  realistic  relationship  can  be  used.  A  great  deal  of  flexibility  can  be  added  to  an 
analysis  code  by  using  a  tabular  method  for  describing  the  constitutive  relationships.  No 
restrictions  are  then  placed  on  how  complex  the  relationships  can  be.  A  tabular  method  for 
describing  the  constitutive  relationships  is  used  in  GROWFLOW. 

The  tabular  method  used  in  GROWFLOW  consists  of  representing  the  constitutive 
relationship  data  with  a  piece-wise  linear  fit.  Figure  5  shows  how  the  piece-wise  linear  fit  is 
made  to  the  constitutive  relationship  data.  Any  desired  point  on  the  curve  is  found  by  linearly 
interpolating  between  two  points  on  the  piece-wise  linear  fit  curve. 
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SECTION  VI 


SUMMARY  AND  CONCLUSIONS 


A  general  purpose  analysis  code  has  been  developed  to  simulate  variably-saturated 
groundwater  flow  in  a  porous  media.  The  code  is  called  GROWFLOW  and  is  based  on  the  SPH 
analysis  method.  An  improved  SPH  method  called  “Gradient  SPH”  was  developed  for  this 
effort.  Simulation  around  boundaries  and  discontinuities  is  dramatically  improved  with  this  new 
method.  Gradient  SPH  is  a  Lagrangian  method.  The  governing  partial  differential  equations 
were  cast  in  a  Lagrangian  frame  of  reference  for  this  analysis  method.  The  governing  differential 
equations  were  then  cast  into  the  Gradient  SPH  formalism.  An  explicit  time  integration  method 
is  used  in  the  GROWFLOW  code,  and  is  therefore  conditionally  stable.  A  formulation  for 
calculating  the  stability  limit  (largest  allowable  time  step)  for  the  explicit  time  integration 
method  is  presented.  GROWFLOW  uses  a  highly  flexible  method  for  controlling  flow  at 
external  boundaries  and  around  internal  boundaries.  Flow  control  panels  are  employed,  which 
are  possible  due  to  the  Lagrangian  nature  of  the  formulation.  GROWFLOW  has  two  types  of 
source  region  models  available:  constant  head  and  constant  velocity.  A  special  treatment  for 
these  regions  is  presented.  A  method  for  representing  the  constitutive  relationships  with  tabular 
data  is  presented. 
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Explicit  Time  Integration  Loop 


Figure  1.  Solution  Flow  Diagram. 
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Interaction  of  Integration  Points 


Figure  2.  Interaction  of  Integration  Points. 


Ditierentiai  Volume 


Figure  3.  Differential  Volume. 
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Operation  of  Source  Region 


Figure  4.  Operation  of  Source  Region. 
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Figure  5.  Piece-Wise  Linear  Fit  of  Constitutive  Relationships. 
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